#include <stdint.h>
#include <limits.h>
#include <math.h>

#include "myexp.h"
#include "myexp_impl.h"

#if defined(MY_EXP_V1)
#define MYEXP exp_v1
#elif defined(MY_EXP_V2)
#define MYEXP exp_v2
#else
#error must define MY_EXP_V1 or MY_EXP_V2
#endif

#include "myexp.inl"

static double specialcase(double x)
{
    uint32_t abstop = top12(x) & 0x7ff;

    if (unlikely(abstop >= 10 + DBL_EXP_BASE)) // x <=-1024(2^10) or x >= 1024.0
    {
        if (isnan(x))
            return NAN; // return nan
        if (x > 0)
            return __builtin_inff();
        else
            return 0.0;
    }

    double r, r2, r4;
    double scale = 1.0;
    double z = x * INV_LN2; // z = x / ln(2)
    union
    {
        double f;
        int32_t i[2];
    } zi = {z + 0x1.8p45}; // zi.i[0] = (round)(z *128)

    int k = (zi.i[0] >> 7); // k= floor(zi.i[0]/128.0)
    int idx = zi.i[0] - k * 128;

    // here, -1477 <= k <= 1477,0 <= idx < 128
    if (k > 1023)
    {
        scale = 0x1p1023;
        k -= 1023;
    }
    else if (k < -1022)
    {
        scale = 0x1p-1022;
        k -= (-1022);
    }

    //-1477-(-1022)=-455,1477-(1023) =454
    // so -1022 <= k<= 1023

    // z0是z的近似值，z0具有这样的形式，z0 = k + idx/128.0
    double z0 = zi.f - 0x1.8p45;
#ifdef MY_EXP_V2
    r = x - (z0 * Ln2hiN) - (z0 * Ln2loN); // 将ln(2)拆分为２个数，减少误差, r= x - (i+k/128)*ln(2)
#else
    r = x - (z0 * LN2);
#endif

    r2 = r * r;
    r4 = r2 * r2;
    double tail = r2 * (C2 + r * C3) + r4 * (C4 + r * C5) + r;

#ifdef MY_EXP_V1
    double body = get_double(exp_data_52[idx], k);
    double res = (body * tail + body);
#else
    uint64_t n1 = (exp_data_64[idx] >> 12);
    uint64_t n2 = (exp_data_64[idx] & 4095);
    double body0 = get_double(n1, k);
    double body1 = div_2p64(n2) * pow2(k);
    double res = (body0 * tail) + body1 + body0;
#endif
    return res * scale;
}

/* 求exp(x)的核心算法，
 基本原理，
  令 e^x　=2^z            (1)
  两边取对数得 x = z*ln(2) (2）
  =>　       z = x/ln(2)
   则　e^x = 2^z= 2^(x/ln(2))
   我们取z的近似值为z0, z0具有这样的格式, z=k + i/128,这里i,k是整数,0 <= k < 128．
   z0满足 z-1/256 <= z0 <= z+1/256, 即abs(z-z0)<=1/256
   z0的计算过程为：
    1. t = round(z*128.00), 将z放大128倍，并四舍五入整，这里，round(z*128.00)是利用浮点数的特点实现的
    2. ki = t >> 7, t>>7等价于t除以128向负无穷取整，当t是非负整数时，t除以128向0取整等于(t>>7),
      而当t就一个负数时，t除以128向0取整并不等于(t>>7), 如-129除以128向0取整=-1,而 -129>>7=-2
    3. idx= t - ki *128,
    4. z0 = ki + idx/128,下面的代码，对z0的计算是利用浮点数的特点实现的
*/
double MYEXP(double x)
{
    uint64_t abs_u64 = as_uint64(x) & LLONG_MAX;

    if (unlikely(abs_u64 > LN2_X_1022_AS_UINT64)) // x>1022*ln(2) or x < -1022*ln(2)
        return specialcase(x);

    double z = x * INV_LN2; // z = x / ln(2)
    union
    {
        double f;
        int32_t i[2];
    } zi = {z + 0x1.8p45}; // zi.i[0] = (round)(z *128)

    int k = (zi.i[0] >> 7);
    int idx = zi.i[0] - k * 128;

    // z0是z的近似值，z0具有这样的形式，z0 = k + idx/128.0,
    // here, -1022 <= k <=1022, 0 <= idx < 128
    double z0 = zi.f - 0x1.8p45;
#ifdef MY_EXP_V2
    double r = x - (z0 * Ln2hiN) - (z0 * Ln2loN); // 将ln(2)拆分为２个数，减少误差, r= x - (i+k/128)*ln(2)
#else
    double r = x - (z0 * LN2);
#endif

    double r2 = r * r;
    double r4 = r2 * r2;
    double tail = r4 * (C4 + r * C5) + r2 * (C2 + r * C3) +  r;

#ifdef MY_EXP_V1
    double body = get_double(exp_data_52[idx], k);
    double res = (body * tail + body);
#else
    uint64_t n1 = (exp_data_64[idx] >> 12);
    uint64_t n2 = (exp_data_64[idx] & 4095);
    double body0 = get_double(n1, k);
    double body1 = div_2p64(n2) * pow2(k);
    double res = (body0 * tail) + body1 + body0;
#endif
    return res;
}
